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An orthogonalization algorithm for producing the pseudo- 
inverse of a matrix is described, and a FORTRAN program which 
realizes the algorithm is given in detail. 
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To every matrix A there corresponds a unique matrix A + with the following 
properties: 


AA + A = 

A 

(1) 

A + AA + = 

A + 

(2) 

(a + a) t = 

A + A 

(3) 

(m*) t = 

AA + 

(4) 


Penrose [l], one of the originators of this concept, called A + the generalized 
inverse of A, and equations (1) through (4) are often called Penrose's Lemmas. 
Recent usage applies generalized inverse to any matrix satisfying (1), (1) and 
(2), or (1), (2), and (3), referring to the unique A + as the pseudo-inverse of A. 
Other definitions of A + have been given (e.g. Albert [2] , Ben-Israel [3]) but the 
most common is that given above. 

For simplicity's sake, the rest of this paper considers only real matrices, 
although most results hold for complex matrices as well. The pseudo- inverse 
provides a way to handle the ubiquitous matrix-vector equation 


Ax - y . (5) 

If A is square and non- singular, A + is A -1 and the vector A + y solves the 
equation. The particular advantage of the pseudo- inverse appears when A is 
singular or non-square, since A + y then is the minimal vector for this equation; 
that is, if M is the set of all vectors x 0 such that 

||Ax 0 -y|| < U Ax - y |l (6 ) 


for all x , then A + y e M and 

i,A+yi! = ^ji Ax °i|y ■ 

Here we use the standard Euclidean norm. 
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A theorem which dates back to the time of Gauss (Newhouse [4] ) states, in 
effect, that if x Q e M , then x 0 is a solution of 


A T Ax - A T y . 


This type of system, often called a set of normal equations, is found repeatedly 
in least squares problems. (See, e.g., Rao [5]). Since A + y e M, the application 
of A + in these circumstances is evident. 

The same theorem also states that if x 0 e M , then x 0 is a projection of y 
onto the column space of A. Newhouse later gives a theorem which proves con- 
dition (7), that A + y is the "shortest" of these projections, giving rise to 
Greville's assertion [6] that A + y is the best solution to equation (5) in the least 
squares sense. 

Naturally, the theoretical existence of such a useful mathematical object 
makes a method for its computation very desirable. Most of the methods sug- 
gested, however, require that the product A T A be formed and that Gaussian 
elimination (or one of its variants such as pivotal condensation or sweep out) 
be performed on it. Should we be faced with an ill-conditioned matrix, it is 
entirely possible that numerical difficulties will prevent any significant com- 
putation using such methods. For example, consider the matrix 


1 1 1 1 

5 6 7 8 

1111 

6 7 8 9 


7 8 9 10 

I 1 1 J_ 

8 9 10 11 


The Hilbert matrix is notoriously ill-conditioned with respect to Gaussian 
elimination. The upper left-hand 4x4 corner of it has a condition number 
X-maxAmin & iven by Marcus (Ref. [7] ) as 15,514, so our 4x4 segment of it 
would certainly be suspect. Fox (Ref. [8]) shows that our suspicions are 
justified, giving the Gauss elimination process forH , in which the steady de- 
crease in magnitude of the pivots leads to very unreliable quantities. More to 
our point, he demonstrates that Gauss elimination fails completely when applied 
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to H T H. We should realize that a bad but workable problem can become path- 
ologically unmanageable if such a product is formed, and, as a general rule, 
avoid such approaches. 

The method of Rust, Burrus, and Schneeberger (Ref. [9] ) was used to 
compute the pseudo- inverse because it does conform to this general rule. 
Briefly, it can be characterized as follows: if the m x n matrix A is in the form 
[R | S] , where the k linearly independent columns form the submatrix R and the 
linearly dependent columns form the submatrix S , make up the n x n identity 
matrix and write, symbolically, 


R | S 

h o 


Then perform the Gram- Schmidt (G.S.) orthogonalization process on [R t S] , 
and apply these elementary column operations to the lower submatrix to get 


Q 1 0 


Z -U 

_° X n- k _ 


Next, perform the G.S. process on the submatrix 



to produce 


Q 1 0 
Z -UP . 

0 P 


form the matrix 


Q t 

(UP) t ZQ t 
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and, finally, 



-UP 

P 


Q t 

(UP) t ZQ t 


A complete derivation is given in Ref. [9] and a few auxiliary notes are given 
in Appendix A. 

Of course, not every matrix will be in the convenient [R| S] form, but if we 
can determine which columns of A are dependent we can certainly permute 
columns to produce it; then [R|s] + is found and by the authority of Theorems I 
and n, Appendix A, the rows of [R| S] + are likewise permuted to get A + . Since 
the G.S. process not only orthogonalizes the independent columns of A but also 
makes the dependent ones zero, we can use it to find the dependent columns. 


Now we have a straightforward way to proceed: 

(1) Use G.S. to find the dependent columns. 

(2) Permute to get [R|S] . 

(3) Use G.S. to find [R |S] + . 

(4) Permute to get A + . 

The reader will have noticed that the G.S. process is used in step (1) and again 
in step (3). We could save some computation time if we combined the two steps 
and performed the G.S. process only once. A closer examination of the process 
reveals that we can, under certain conditions, make this combination. 

Our program uses a modified Gram-Schmidt process which is more accurate 
than the classic textbook version. A recursive algorithm describing our version 
is: 

(1) Orthogonalize c. , the next column of A: 


b. - c - 

j j 



b ' 
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(2) Is bj % 0? If so, zero it out and go to step (1). If not, do step (3). 

(3) Re-orthogonalize b j : 



and go to step (1). 

The initial condition is b x ' = c x . After we run out of columns, we normalize 
each one and we have an orthonormal matrix A x . 

If we want to duplicate these elementary column operations on another 
matrix D, we could save the numbers 



and (bj'-b/) 1 / 2 


and then go through the algorithm again, this time letting be the columns of 
D. More precisely, we might save these numbers in an n x n matrix S , defined 
as 


c . • b 

j i 


V 

ii 

cr 1 j 

i 

• b.' 

(j > i) , 

s.. 

i 

_ 

1 

■v 

(j > i) , 


b.' 

• b.' 


i 

1 


s ,i 

•n" 

ii 

i'-B/r 

( 1 < j in) 


As an example, let A be 

”l 1 3 6 

2 2 6 7 

3 3 9 8 
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Using eight-digit arithmetic and rounding the final answers to three digits, 
we have 


213 x 10' 


Once we have done the G.S. process on A, we have done it for all column permuta- 
tions of A which do not disturb the relative order of the independent columns. If 
P is a permutation matrix such that AP n [R| S] , where R is the matrix of inde- 
pendent columns of A in their original relative order, the orthonormal matrix 
[Q| 0] produced by the G.S. process on [R|S] will be AP . In our example, 
suppose 


Then 


[Qlol 


-.436 


We can also produce a new S matrix (F) by permutations. Referring back to 
the definition of S , one can see that a particular column c has its initial 


6 


orthogonalization coefficients (c. . b. )/(b L . b. ) on the j th row and below the 
diagonal, and its secondary coefficients (b. . b')/(b' . b'jfallonthe jth column 
and above the diagonal. Once c. is converted into b . ' , all initial coefficients 
having b/ as a factor fall on the jth column below the diagonal, and all such 
secondary coefficients fall on the jth row above the diagonal. When b/ is 
normalized, its ’length" falls on the jth diagonal element. Moving Cj to a new 
position therefore means that we must move the jth row and column of S to 
corresponding positions, or 


F = P T SP 


In our example, 


F 


3.74 

.213 xlO -7 

0 

0 

3. 14 

3.27 

0 

0 

1.00 

0 

0 

0 

3.00 

0 

0 

0 


If we go through the algorithm with c taken as the columns of a matrix D and 

9 9 9 i I 9 9 

the numbers (c. . b i )/(b i . b.)and(b. . b i )/(b i . betaken as f jjL and f.. re- 
spectively, we have applied the elementary column operations of the G.S. process 
on [R| S] toD. 

Now we have the desired result: once the G.S. process on A is complete, it 
is not necessary to do it again on [r| S] to derive its effects; merely execute the 
indicated permutations on Ai and S and we have all the necessary matrices. 

Using this result, the procedure (1) through (4) on page 4 can be rewritten: 

(1) Use G.S. on A; save the G.S. coefficients in S and save A x . Note which 
columns are dependent. 

(2) Permute A x to get [Q| Ol ; permute S to getF. 

(3) Use the entries of F to operate on 
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producing 



(4) Proceed as usual to find [R| S] + . 

(5) Permute to get A + . 

The program whose flow chart and FORTRAN listing appear in Appendices 
B and C has been checked with a variety of matrices on the IBM 7094 and appears 
to run properly. Two particular cautions might be extended, however: first, 
one will note that a decision on the dependency of any column is made by com- 
paring the "length" of the generated orthogonal column with the "length” of the 
original column. If the check number (b. . b. )/(c. . cp is smaller than a certain 
tolerance, the column b f is made zero. When the check number is very close to 
the tolerance, any decision made will not be a good one and the resulting per- 
turbations can become serious; for example, the Hilbert matrix gives poor 
results for this very reason. One might vary the tolerance to suit special cases. 

Second, although this program finds the inverse if it exists, there are 
routines in general use which get better inverses. For example, the SHARE 
routine MATINV was tested against this program on a sequence of Pei matrices 
(Ref. [10], Ref. [11]) and consistently got one more accurate digit in the worst 
cases. The difference is not great but the prospective user should realize that 
it exists. 

Finally, an experienced programmer will see that the FORTRAN realization 
in Appendix C is not in optimal form. A more streamlined, double-precision 
version is being prepared for the IBM 360 as of this writing. The author would 
appreciate hearing of mistakes in, or improvements upon, the original. 
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/ APPENDIX A 

✓ 

I . 

(Supplementary notes for Ref. [9] ) 


Theorem I 

If P is a permutation matrix (possibly a product of elementary permutation 
matrices) and A + is the pseudo-inverse of A, then 


(AP) + = P T A + . 


Proof: 

We need only verify that Penrose's Lemmas hold. Noting thatPP 1 

= pTp 

- I, we have 


(a) 

(AP) ( P T A + ) (AP) 

= AP 

(b) 

P t A + (AP)P t A + = 

P t A + 

(c) 

[(AP) ( P T A + )J T = 

(AA + ) T = AA + = (AP)(P t A + ) 

(d) 

[( P T A + )(AP)] T = 

[ p t (A + a) p] t 


= P T (A + A) T P 

= P t (A + A) p 


= ( P T A + ) (AP) . 


Theorem n 


If P is a permutation matrix and the operation AP effects a column permuta- 
tion of A, then P T A effects that same permutation on the rows of A. 

Proof : Suppose one of the effects of AP is to change column i to the j th place. 
ThenP.. = 1, P. T . = 1, and P T A changes row i to the jth place. 

We use this result to get A + from a row permutation of [r| S] + — that same 
permutation of columns which transformed A into [R| S] . 

The paper states (p. 383, right column) that the G.S. process turns a de- 
pendent vector into the zero vector. One might check this statement by referring 
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to Hoffmann and Kunze, p. 230, Theorem III (Ref. [12]). If a k+1 is a linear 
combination of a t , • • • , a k then it is a linear combination of q lf • • • , q k since 
the vectors q i span the space of the vectors Furthermore, by the above- 
mentioned theorem, 


a 


k + 1 




di 


and c k+1 n 0. 

On p. 384, left column we are to note that I n _ k remains unchanged. Suppose 
we are operating on column k + p (p>0) of the matrix [R| S]. We have 


k+p-1 


' k +p 


a 


k +p 


£ ( a k" P <ii) q, 

i = l 


k k +p- 1 

= 3k+ p ■ 2Z ^ ak " pq ^ " L ( akH+ p q ^ qi * 

i= 1 i=k + l 

But each q.,k+l<i<k + p-l, has been zeroed out already, since they came 
from vectors dependent upon a a k , so the above is 


c 


k +p 



Similar column operations on the identity matrix then use only the first k columns, 
whose lower n - k entries are all zero and cannot contribute to any modification 
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APPENDIX C 

“uTPIINbi on m loVicn » factort lo» rol^5T:rWf5t^^ 

lAlCIQ? IQ) ?ORI 0 ( 10,10) , PROD (10,10) , PRQD1 ( 1 Q , 10 ) , L __ 

2PRGD2 (10*10) »ATR(1Q»1G) *AINV (10*10) 

smm&Uth- uptr ( 10 , 10 ) . 

DIMENSION AINV1 ( 10 » 10 ) 

D A T A PR A IN1/5H A INV1/ — 1 

DATA PRA»PRFACT »PR52/1HA*4HFACT *2HS2/ . I 

; DA T A PRA1/2HA1/ ; ; — ! : - ! , 

DATA PRAINV/4HAINV/ 

Dm I A. PkuR I u »PRPR « ERBxLL « PR P A/ .PR , RTk/4nutiC I Q, 4 HPRQD . 5HPRQD1 .5 HP Ron? f ; . 

13HATR/ 

DmTA PRLiPTk/ 4KuPTR / . ;j 

134 READ (5,1) NROWS »NCOLS | 

2 FORMAT!? 15) _ : ! 

DO 112 I = 1 , NROWS / _ . ;i 

2-12 READ! 5 .2 ) .I.AII,.J).J = l.NCOLS) _li 

2 FORMAT (6E12. 8) 

T 01 = t 1 n.*n. 5dULZ2 1 *±2 __J 

DO 110 I = 1 -.NROWS ■ ' 

DO 110 J = l.NCOLS ' - j 

110 OR 1 6 ( I » J ) = A ( I » J ) , ’ 

DO 100 1=1 ,NCOLS - 

I HOLD ( I ) * 0 

JHOLD ( I ) = 0 : ; 

DO 102 J = 1 » NCOLS “ • ~ \ 

S2 < I § J ) = 0. 

PROD1 ( I » J ) = 0. 

PROD2 ( I , J) = 0. | i 

ATR ( F7J) = 0. ' 

FACTOR ( I , J ) = 0, 

I 02 All I *J) = 0. 

loo a i ( i » i ) = l. .4 ; - 

J H0LD(1) _ l ‘ ; 

kk = i ; - : 

_ ■ T ! : 

11=0 , . ! ' 


152 NLESS1 = N 

N = N + 1 

CHECK = DOT (A«N.N. NROWS) _____ 

DO 101 I = 1 »NLESS1 

F ACTOR (N,l) = DOT I A.Ni 14 NROWS ) /DOT ( A , I « I . NROWS ) 

DO 101 J = 1, NROWS 

HU A ( J » N ) = A(J,N) - F AC IQR ( N , 1 ) * A ( J , I ) - 

CHECK = DOT ( A»N »N »NROWS ) /CHECK 

, I F ( CHECK - TOL) 150,150,151 ~ 

150 DO 103 J = 1, NROWS 
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JJ = JJ + 1 

IHOUD(JJ) = N ' 


GO TO 155 

151 

DO 104 I = ltNLESSl ! 

F ACTOR ( I * N ) = DOT(A»N*I *NROWS)/DOT(A*I»I» NROWS ) , 
DO 1C4 J = 1, NROWS 

104 

A ( J » N ) = A ( J »N ) - FACTOR! I »N)*A(J, 1 ) i 

KK = KK + 1 ..... . ] 


JHOLD ( KK ) = N 

155 

I F ( N - NCOLS) 152*153*153 1 

| 

L 

DO 105 J = 1 *NCOLS 

FACTOR(J»JJ = SORT! DOT ( A*J*J*NROWS) ) , 

I F ( FACTOR ( J *J ) « EQ • 0.) GO TO 105 ! 


DO 106 K = 1 * NROWS 

106 A ( K * J ) = A(K,J ) /FACT OR ( J aJJ 

ic 5 CONTINUE 

LAI I WRITE (FACT OR »NC01 S.roi S.PRFACT ) 

CALL WRITE('a,NROWS*NCOLS,PRA ) 

IF (KK.EQ. NCOLS I GO TO 1 5 A 

DO 120 I = 1 * KK 

I SU B = J HOL D! I) 

DO 120 J = 1 * NCOLS 

A(J*I ) = A(J*ISUb) 

120 S2 ( J * I ) = F ACTOR ( J » I SUB ) 

DO 171 T = 1 .KK 

I SUB = JHOLD ( I ) 

D O 1 2 1 J. = 1 » NCOI .S 

121 FACTOR ( I * J ) = S2(ISUb»J) 

K K = KK + 1 - 

DO 125 I = KK * NCOLS 

LL =—.I. 1 . . ±-.3 

ISUB2 = I-HOLD! I I ) 

DO 12 5 J = I » Nf 01 S 

125 FACTOR! I * J ) - S2!ISUb2,J) 

- DO 162 I a K K * NCOLS 

DO 162 J = 1 * NROWS 

.162 A ( .,1 » I ) a 0. 

KK = KK - 1 

166 CALL. W R I T F ( A > NROWS * NCOLl-»-PRAJ 

CALL WR I TE ( FACTOR , NCOLS* NCOLS *PRF ACT) 

A L AI I W R 1 T E ( A 1 . NCOLS *NLU_LSaP R A1 .) 

DO 502 I = 2, NCOLS 

. I LESS 1 = I -1 

DO 500 J = 1*1 LESS 1 
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90 500 r.« UNROWS l 


a 

1 

!•* I* rdHsaV V ■ A ■ mw ■KLJHV* P V m Cm ■ W Haim M m ■ ■ m mPm ^ — 1 

DO 501 1*1 »KK 
DO 501 J * 1 , NROWS . 

501 A 1 ( J , I ) = A1 ( J * I) /FACTOR ( 1,1 ) 1 

CALL . wft I TE ( Al » NCOLS »NCOLS » PRA1 ) ■ • -v 

CALL TRANSP ( A » NROWS » KK , ATR ) i 

CALL WRITE(ATR,NCOLS, NROWS, PR AIR ) j 

. I F ( KK • EQ , . NCOLS ) GO TO 157 

1ST a KK + 1 1 

IF(KK.EQ. (NCOLS-1 ) ) GO TO 158 ' ! 

.‘ v N a Kk + 1 x • 

159 i NLESS1 * N 

_ N = N +1 _ . 

• DO 107 1 t 1ST »NLESS1 

FACTOR JN. I » = DOT ( A1*.N « 1 • NCOLS 1 /DOT 1 A 1 • I 1 1 »NCOLS l 

DO 107 J a 1, NCOLS 

4,07 AT ( J f N t a AlfJ.N) - E ACTOR ( N «_L 1 * A1LJ f I) 

DO 108 I = 1ST »NLESS1 

FACTOR < I .N) = DOT i Al . N . I . NCOI S \ /DOT l f 1 . I f I ,NCOI S V 

DO 108 J a 1, NCOLS 

1 OR A 1 ( J j N > a Al t J.N) - FACTOR ( I ,N ) * AI(J.I) | 

IFIN.LT. NCOLS) GO TO 159 
.CALL wRtTtTAl .NCOI S,NCOI S.PRA1 ) 

158 DO 128 1 a 1ST »NCOLS 

U^. FACTOR ( t » I ) g SORT ( DOT ( Al »t . l.tNCOI S ) ) ~ 


DO 128 J * 1 ♦ NCOLS 


T2r"""7rnT»Il~='ATTG7I T/F 5CWn7TT^ - ~ 

i I , C A L L W Kl T . LtA UNCOL^M-CQLS tPR A l) , 

CALL MATMPY < Al »KK-»KK » ATR »NR0WS*PR0D1 ) 

■ » ■ »» CALL WR 1 Tfc( PROD1 > KK > NROWS , PRPR1 j 

LIM = NCOLS - KK 

— DO - 122 I - 1 > LI M 

ISUB1 = KK + I 

f- 00 122 J - » 1 > NCOLS : 

122 UPTR ( I » J ) = -Al(J,ISUbl) 

I CALL WRITE ( OPTR »- L I M » NCOLS , PR U -P- T -R K- 

CALL MATMPY ( UPTR » L IM » KK »P ROD 1 » NROWS, PROD 2 ) 

C All URTTF<PROn?.l T M . NROWS . PR PR ? 1 

DO 123 I = 1*LIM 

■- 1SUB3 = KK ± t 

DO 123 J = 1* NROWS 






CALL WRITE(ATR»NCOLS, NROWS* PRATR ) 

157 CALL MATMRY( A1 » NCOLS , NCOLS » ATR » NROWS » A I Nw ). 

CALL WRITE! A1NV»NCOLS»NROWS»PRAINV ) 

CALL WRITE(ORIG*NROWS > NCOLS » PROP I G ) 

DO 126 I = 1 » KK 

ISUBA = JHOLD(I) 

DO 126 J * 1 * NROWS 

126 A I NV1 ( I SUB4 i J ) = A I N V ( I » J ) : 

IF (KK .EQ. NCOLS) 00 TO 160 

K K a KK + 1 

11=0 

00 1 ?7 T = KK .NCOl S 

II = II + 1 

I SUB 5 = 1 HOI Dill) - 

DO 127 J = 1 *NROWS 

127 A I NV 1 l I SUBS » J ) = A1NV(I»J) 

KK = KK - 1 

160 CALI. WRITL(AINV1 *NC0LS»NR0WS»PRAIN1) 

CALL MATMPY (OR I G ( ,NROWS> NCOLS* AINVl , NROWS, PROD ) 

CALL Wk I TE( PROD *NROWS* NROWS *PRPR) 

CALL MATMPY ( PROD, NROWS* NROWS > OR I G * NCOLS » PR0D2 ) 

CALL WRITE(PR0D2*NR0WS» NCOLS * PRPR2 ) . 

CALL MATMPY ( AINVl , NCOLS * NROWS* PROD * NROWS, PR0D2 ) 

CALL WR I TE ( PR0D2* NCOLS* NROWS >PRPR2 ) 

CALL MATMPY (AINVl , NCOLS ♦ NROWS * OR I G ♦NCOLS, PROD ) 

CALL WRITE! PROD *NCOLS*NCOLS*P.RPR) 

GO TO 154 

EM) __ 


FUNCTION DOT ( A, J.,K*NROWS ) 

DIMENSION AMn.10 1 

DOT =0.. 

DO 1 0 I a 1 , NROWS 

10 DOT = DOT + A(I.J) * A ( I » K ) 

RETURN 

END 


I .SUBROUTINE WR IT LI X , NRO WS « NCOLS * NAM E ) 

DIMENSION X(10.10) 

; WRITE ( 6 « 3 ) NAME 

3 FORMAT ( 1H0////7H MATRIX»1A7) 

DO 21 I = 1 * NROWS 

21 WRI TE ( 6 *4 I (X(I»J)»J = 1, NCOLS) 
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4 FORMAT ( 1H0»8E16*8) 

RETURN 

END 


SUBROUTINE MAT-MPY ( A « NRA «NCA « R « NCB . PROD) j 

DIMENSION A(10»10)»B(10'*10)» PROD ( 10*10) | 

DO 600 -J = 1»NR,A — J 

DO 600 J = 1 »NCB | 

PROD ( I ». J > =— 0 . I 

DO 600 K = 1 »NCA | 

£00 PROD ( I » J 1 . 5 - PR-QD< !-»■))-. + A(T»K) £ R ( K » ..I ) j 

RETURN | 

EJSIO : 1 



SUBROUTINE TRANSP1 A * NRA »NCA , ATR ) 
DIMENSION AT 10*10 L*AT_RJ 10*10) 


DO 601 I = 1 *NRA 
_DO 601 J = 1 ,NCA 

601 

ATR ( J » I ) = A ( 1 * J ) 
RETURN 


END 
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